The Dynamics of Sustained Reentry in a Loop Model with Discrete Gap Junction 
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Dynamics of reentry are studied in a one dimensional loop of model cardiac cells with discrete 
intercellular gap junction resistance (R). Each cell is represented by a continuous cable with ionic 
current given by a modified Beeler-Reuter formulation. For R below a limiting value, propagation is 
found to change from period- 1 to quasi-periodic (QP) at a critical loop length (L cr it) that decreases 
with R. Quasi-periodic reentry exists from Lcrit to a minimum length (Lmin) that is also shortening 
with R. The decrease of L cr i t (R) is not a simple scaling, but the bifurcation can still be predicted 
from the slope of the restitution curve giving the duration of the action potential as a function of 
the diastolic interval. However, the shape of the restitution curve changes with R. 

PACS numbers: 87.19Hh, 05.45-a 



I. INTRODUCTION 



Self-sustained propagation of electrical activity around 
a one-dimensional (1-D) loop of cardiac tissue is the 
simplest model of reentry, the mechanism by which a 
propagating activation front maintains itself by travelling 
around a functional or anatomical obstacle. Reentry has 
been much studied because it was demonstrated to be an 
important mechanism of cardiac arrhythmia. [J 0, d, H| 
For the 1-D loop, most work has been done assuming 
the membrane to be a continuous and uniform cable 
with constant intracellular axial resistivity. 0, IE 0, H, M, 
EE EL GJ, EL EI E| For different models represent- 
ing the ionic properties of the membrane, propagation 
was found to change from stable period- 1 propagation to 
quasi-periodic reentry when the length of the loop was 
reduced below a critical length. The quasi-periodic reen- 
try was characterized by a spatial oscillation of the action 
potential duration as propagation proceeded around the 
loop. Based on numerical simulation, the bifurcation was 
in most cases classified as supercritical, with the ampli- 
tude of the oscillation growing as the length of the loop 
was reduced below the critical length. Quasi-periodic 
reentry was found to exist from the critical length to a 
minimal length below which sustained propagation be- 
came impossbile. In some instances, two different modes 
of quasi-periodic propagations were identified, with dif- 
ferent wavelengths, different interval of existence, and 
sometimes different scenarios of creation. [EBBS EE] 
Different attempts were made to build simplified repre- 
sentations of the dynamics allowing analytical examina- 
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tion of the nature of the bifurcation. @, El, El EE EL El 
One of these approaches, which guides the present inves- 
tigation, relies on an integral-delay model. [E EE It is 
based on the assumption that both the speed of propaga- 
tion and the action potential duration can be expressed as 
functions of the diastolic interval, which measures the re- 
covery time from the end of the previous action potential. 
The model has been successful in reproducing the locus 
of the bifurcation observed by numerical simulations of 
1-D loops with Beeler-Reuter type representation of the 
membrane. It predicts that the bifurcation should occur 
when the diastolic interval in the period- 1 reentry reaches 
the critical value where the slope of the restitution curve 
becomes 1. 

However, cardiac excitable tissue is not a syncytium, 
but rather a mesh of myocytes connected by discrete gap 
junction resistance.: 19] Much work has been done to in- 
vestigate the effect of discrete gap junction resistances 
in a one-dimensional structure, |lJ |2E E3, IHIM 13 
mostly focused on the effect of resistivity on excitability. 
In the discrete case, the resistance no longer acts as a 
scaling factor with regard to space. Because the inter- 
cellular current is reduced as the gap junction resistance 
is increased, the latency of the cell to cell propagation is 
increased until propagation fails at some limiting value 
of the resistance. Besides, upon premature or repetitive 
stimulations, the excitability of the tissue must be more 
recovered for propagation to proceed, which corresponds 
to an increase of the refractory period. 

This paper describes how the bifurcation from period- 
1 to quasi-periodic propagation and the characteristics 
of the quasi-periodic propagation are modified by the in- 
crease of the intercellular resistance in a 1-D loop of dis- 
crete model cardiac cells. This paper is organized as fol- 
lows. In the next section, the model and computational 
method are described. The results of the numerical sim- 
ulation are presented in §3. The bifurcation from stable 
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period-1 reentry is explained in §4. The QP modes of 
reentry are analyzed in §5. The final section is devoted 
to a summary and discussion. 



II. METHODS 

We consider a one-dimensional loop formed by N iden- 
tical cells connected by gap junction resistanced. Each 
cell is modeled as a continuous and uniform cable of ra- 
dius (a) 5 pm, length (L c ) 100 pm and intracellular re- 
sistivity (p) 0.2 KQ-cm lying in an unbounded volume 
conductor of negligible resistivity. The transmembrane 
potential (V i=1,N in mV) of the cells is described by the 
well-known cable equation: 



1 <9 2 ^(a;,i) 
P 



S\C, 



0V*(a;,t) 



dx 2 dt 

x e {0,L C }, it {1,7V} 



(1) 



in which C m is the membrane capacitance (1 pF/cm 2 ), 
S is the surface-to- volume ratio (0.4 pm^ 1 ) and I ion is 
the ionic current ( pF/cm 2 ). The membrane ionic model 
is the same MBR model that was used in ourprevious 
works on continuous 1-D and 2-D rings. @, 0, IEEE El 
In this model, the sodium current is controlled by two 
inactivation gate variables h and j . Each cell is connected 
to its neighbors by a discrete gap junction resistance R 
(Ktt). Continuity of the intracellular current between 
the cells yields the boundary conditions: 



dV l 



Qymod(i,N) + l 



. . x—L c . . 

OX OX 

V*(L C ) - V m ° d ^ +1 (0)=RI hmodihN)+1 (2) 

For simulation, we have modified the numerical 
method that was developed for continuous loops. [6| 
Briefly, for each time step (At = 2 ps), Eq.([T]) becomes 
equivalent to an ordinary differential equation 



\ x=0 — 9 Ii,mod(i, JV) + 1 



d 2 V l {x) 
dx 2 



K V" 1 (x) = g\x) 



(3) 



whose solution can be expressed as the sum of a partic- 
ular solution Vp (x) and of the homogeneous solution 



The purpose of the simulations is to obtain a descrip- 
tion of the regimes of reentry of the function R and 
L = NL C , the length of the loop. During reentry, the suc- 
cessive action potentials (j = l,m) at each node can be 
characterized by their activation times {T ] act ), set at the 
maximum derivative of the upstroke, and their repolar- 
ization times (T 3 epol ), taken at the —50 mV downcross- 
ing in repolatization. The action potential duration (A) 
and the diastolic interval (D) associated to each action 
potential are calculated respectively as A> = Tf 



repol 



-T 3 

act 



and D J = T 3 act - T° r ^ ol . @, i, 0, 1 The propagation of the 
wavefront along the loop generates spatial profiles of A 
and D that typify the reentry. In contrast to a continu- 
ous loop, propagation on a discrete loop can be patterned 
inside each cell but identical across all the cells. We have 
chosen to use only the A and D values of the middle node 
of all cells to characterize the reentry. We label period-1 
(P-l) reentries in which A and D remain constant across 
all the middle nodes, and quasi-periodic reentries where 
A and D oscillate both in time and space. 

For each value of R, an initial L was chosen large 
enough to sustain P-l stable reentry. Reentry was ini- 
tiated by transiently opening the loop and stimulating 
one end. Computation was continued until stable period- 
1 reentry was detected, the stability criteria being less 
than 0.5 ms difference in A and D between all the mid- 
dle nodes for one rotation of the front. Afterward, the 
loop length was gradually reduced by steps of one cell, 
using the final state of the previous L as initial condi- 
tion and removing one cell far from the position of the 
excitation front. When the stability criterion was not ful- 
filled after a minimum of 25 turns, reentry was labelled as 
QP. With this procedure, both L cr n and P cr it, respec- 
tively, the minimum length and minimum period with 
P-l reentry, as well as L m i n , the minimum length for 
sustained reentry, were identified for each value of R. In 
some instances, bistability between P-l and QP reentry 
was investigated by stepwise expanding loops that were 
initially in quasi-periodic regime. One cell was inserted 
in the loop, with initial conditions set at the mean of the 
states of its neighboring cells. Finally, we also searched 
for distinct modes of QP reentry using the method de- 
scribed in [7j, in which the D spatial profile of a QP 
solution for a given L value is compressed by a scaling 
factor and used to construct an initial condition to find 
alternative QP solutions with smaller wavelengths. 



V£(x) = Aie** + Bte-** (4) 

Vp (x) is obtained by solving Eq. §5§ with Neumann 
boundary conditions (dV i /dx\ x — q.l„ — 0) using the lin- 
ear finite element method[13| with a uniform grid ( Ax = 
25pm), i.e. 5 nodes. Cells are then reconnected by 
choosing the coefficient of the homogeneous solutions to 
fulfill the continuity conditions Eq.([2]). For a subset of 
R values, calculations repeated with Ax = 12.5pm and 
At = lps gave the same results. 



III. RESULTS 

Figure [T] A) shows L cr a and L m i n as a function of R. 
Both L cr it and L m j n decrease until they merge at R = 
104 Kfl. From this resistance, QP reentry does not ex- 
ist anymore and P-l reentry remains the only regime of 
sustained propagation. From there, the limiting length 
for P-l reentry increases until sustained propagation be- 
comes impossible at R = 108A29KH. Increasing the 
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resistivity in a continuous loop would also decrease L cr n 
and L min . However, the speed of propagation being pro- 
portional to 1/ y/p in a continuous media, [27| y/pL crit (p) 
and v /pi m i n (p) would remain invariant. To compare the 
continuous and discrete medium, we computed the equiv- 
alent resistivity of the latter as 



NRira 2 
Peqv{R) = pH j = P 



Rna 2 



(5) 



If the two media were equivalent, the ratio 
L(R) \J Peqv(R) I 'L(0)y/p would remain equal to 1. Fig. 
[T]B) shows clearly that the diminution of L cr n and L m i n 
cannot be explained by a simple scaling, as it occurs in 
a continuous medium. 
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FIG. 1: A) L cr u (cm, dashed line), the shortest L with period- 
1 reentry, and L m i„ (cm, continuous line), the mininum L 
with QP reentry, as a function of the gap resistance R(Kfl). 
B) Normalized values of Lcrit and L m in (see text) as a func- 
tion of R. 



For each value of R, we collected the D and A values of 
the P-l solutions for a collection of L values as well as 
those of the first QP solution below L cr n to construct 
the A(D) restitution curve. Each curve was fitted with 
a simple exponential to find D cr i tj th{R), the value where 
the slope of the fitted A(D) — 1, and theoretical value 
P C rit,th = D crit .th + A(D crit<th ). As shown in Figd A), 
Pcrit,th falls very close to the P cr it values found by simu- 
lation. Hence, the mechanism responsible for the transi- 
tion from P-l to QP reentry is the same in the continuous 
and discrete loop, and the increase of P CT it results from 
R transforming the restitution curve. The mechanisms 
responsible for the change of D and A can be identified 
in Fig. [31 which shows the action potentials of the first 
node in the three successive cells for increasing values 
of R. (top to bottom, R = 0, 80 and 103 Kfi). In- 
creasing R prolongs the latency of the action potential 
(left column panels). Because the end of the diastolic 
interval is set at the maximum upstroke derivative, the 
increase of latency is translated as an increase of D. On 
the other hand, neighboring cells exchange current dur- 
ing the repolarization, which compensates for the delay 
of activation and prolongs the action potentials. These 
two effects contribute together to the change of the resti- 
tution curve. Finally, the transition from period-1 to QP 
solution was always found to be supercritical, as in the 
continuous case. 



IV. 



P crlt IN TRANSITION TO QP 
REENTRY. 




Two distinct scenarios can lead to the disappearance 
of P-l reentry. For 108.429.fm > R > 104KQ, sus- 
tained reentry does not exist for L < L crit = L min , 
so that reentry ends abruptly with the disappearance of 
the P-l solution. For R < 104KT2, P-l reentry is re- 
placed by QP reentry that persists from L cr u to L m i n . 
In this section, we consider the second type of transition. 
In the continuous MBR loop, the bifurcation from P-l 
to QP propagation was shown to occur at the critical 
period P crit = D crit + A crit where D crit and A crit are 
the values for which the slope of the restitution curve 
A(D) reaches l.Jg, || P cr it is constant and independent 
of p in a continuous medium. In contrast, Fig. [2] A) 
shows that P cr n increases with R in the discrete loop. 
Both A cr it and D cr u contribute to the change of P cr it 
(Fig[2]B), but the increase of D cr n is more important. 



FIG. 2: A) Continuous line: For each value of the gap re- 
sistance R, the critical cycle length Pcrit (ms) at Lcrit, the 
shortest loop with period-1 reentry. Dashed Line: P C rit,th, 
the critical cycle length computed from the restitution curve 
(see text). B) Value of the action potential duration (A cr it, 
continuous line) and diastolic interval (Dcrit, dashed line) at 

Lcrit - 

Once Pcrit is known, L cr n can be calculated if the 
speed of propagation Q(D cr it) is provided. In the dis- 
crete media, the total time to propagate from one cell to 
another is a composite of propagation time within and 
between the cells. The former decreases with R , while 
the latter, which is equivalent to the latency displayed in 
Fig. [31 increases. The final composite Q(D crit ) is shown 
in Fig. Q] 
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FIG. 3: Action potentials (mV) in the first node of 3 succes- 
sive cells as a function of time (ms) during period-1 reentry 
for, from top to bottom, R = 0,800 and 103 KQ. On the 
left column panels, only the activation is shown, while the 
complete action potentials are displayed in the rigth column 
panels. 





FIG. 5: Characteristics of the mode-0 QP solution at L m i n 
for R = 3KQ (top row panels) and R — 103KQ (bottom row 
panels). Left panels: D, the diastolic interval) as a function 
of position (x/L) for 3 successive turns abutted end-to-end. 
Middle panels: A, the duration of the action potential, as 
a function of D, from the solutions shown in the left panels. 
Rigth Panels: 1/0, the cell-to-cell conduction time, as a func- 
tion of D from the solutions shown in the left panels. Only 
the data of the first node of each cell were used to construct 
these plots. 
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FIG. 4: Intercellular activation speed (calculated between the 
first node of successive cells, cm/ sec) at L cr u as a funtion of 
the gap resistance R. 



V. QP REENTRY 



The characteristics of the QP reentry in the continuous 
MBR loop have been extensively discussed in previous 
papers. 0] Two modes of QP were identified, character- 
ized by D oscillations with different spatial wavelengths 
(A). The first mode, referred to as mode-0, exists from 
L cr it to L m i n . Its A, close to two turns of the loop at 
L cr n, diminishes as the loop is shortened, but always re- 
mains longer than L. The second, referred as mode-1, 
exist only over a subset of the [L m i n , L cr i t ] interval with 
a A always less than L. These two types of QP solu- 
tions were found for all values of R < lOAKtt where QP 
solutions exist. 



A. mode-0 QP reentry 

We first consider the mode-0 solutions that exist over 
the whole [L m i n , L cr i t ] interval. Fig. [5]presents the char- 
acteristics of the mode-0 solutions at L m i„ for two val- 
ues of R (top panels , R — 3K£l, L = 7.65cm , bot- 
tom R = 103Kfl , L = 1.04cm ). The leftmost pan- 
els show the spatial oscillation of D by plotting succes- 
sive turns end to end. The first obvious difference is the 
range of D values covered by the two solutions. The 
left panel of Fig. [5] shows D m i n and D maxi the mini- 
mum and maximum value of D for the mode-0 solution 
taken L m i n as a function of R. It is well known that, 
in a discrete medium, the minimum excitability needed 
to sustain propagation increases as a function of R un- 
til a limiting value beyond which propagation is blocked 
even in a medium at rest.[22| In the MBR model, the 
excitability can be measured by the product hj of the in- 
activation gates of the sodium current. The right panel 
of Fig. [S]shows hj(D m i n ) and hj(D max ), the excitability 
of the action potentials produced respectively at D m i n 
and D max for the mode-0 solutions at L m i n . As R in- 
creases, the minimal excitability allowing propagation 
becomes higher, which requires an increase of D m i n (R). 
At R = lOiKtt, hj{D min ) = hj(D crit ), QP propaga- 
tion disappears and only P-l reentry remains. On the 
other hand, the curve hj(D max ) rather reflects the in- 
activation of the socium current occurring during the 
latency preceding the upstroke of the longer action po- 
tential. At R = 104Kfl, the limit for QP propagation, 
hj(D max ) > hj(D mm ) = hj(D crit ), which indicates that 
P-l propagation is still possible if R is increased. How- 



ever, the difference is smaii, such that the range of R 
values over which P-l reentry can still occur is limited, 
as it is seen in FigfT] 



- *j 



FIG. 6: Left Panel: D m i n (continuous line) and D max (dashed 
line) respectively the minimum and maximum diastlolic inter- 
val of the mode-0 QP solutions at L m i n as a function of R. 
Rigth Panel: hj, product of the sodium current inactivation 
gates taken at D m in (continuous line) and D max (dashed line) 
for the mode-0 QP solutions at L m in- 



The middle column panels of in Figj5]display the A(D) 
relation obtained from each QP mode-0 solution. Each 
curve has two branches, the lower and upper branch com- 
ing respectively from increasing and decreasing portion 
of the D spatial profile. Such a dual structure has been 
observed in the continuous loop and was explained either 
by the influence of neighbors on the repolarization [181 ] or 
by short term memory [15j | . The separation between the 
branches is enhanced by the increase of R. Finally, the 
right column panels of Figj5]show 1/9 vs D, the disper- 
sion relation of the conduction time. For R — 2>KVi (top 
right panel), the dispersion relation appears as a single 
value function, similar to what is seen in the MBR con- 
tinuous loop. For R = (bottom right panel), the 
dispersion relation has two branches, as the A(D) curve. 
The lower branch is associated with the decreasing por- 
tion of the D spatial profile. 




FIG. 7: Mode-0 (top panel) and mode-1 (bottom panel) QP 
solutions for L =1.80cm and R =50K£l. The plots show D, 
the diasolic interval, as a function of position (x/L) for three 
turns abutted end-to-end. 



A(n) 



2L 



C 



2n + l (2n + 1) E 



(6) 



where n is the order of the mode and C is a small positive 
constant. As seen Figf?) A(0)/A(1) is indeed close to 3. 
However, in the MBR continuous loop, only the first 2 
modes (i.e. and 1) were observed. This was explained 
by the effect of resistive coupling between neighbors that 
limits the spatial gradient of voltage and forbids the ap- 
pearance of higher modes. Q Theoretically, 



A(2) 
A(0) 



A(2) 
A(l) 



(7) 



To look for mode-2 solutions for different R and L val- 
ues, we have compressed the D profiles of the mode-0 
and mode-1 solutions up respectively to a factor 6 and 
2 to build different initial conditions. This procedure 
was successful in obtaining mode-1 solutions from mode- 
solutions, but higher modes of propagation were never 
produced. For all scaling factors, propagation was found 
to stabilize either to mode-0 or mode-1. 



B. Higher QP Modes 

FigfT] shows an example of mode-0 and mode-1 solu- 
tion for R = 50Kfl and L = 1.8 cm, in the middle of the 
[Lmim L cr it] =[1.61cm, 1.99cm] interval for this value of 
R. Mode-1 solutions were found for all value of R with 
QP propagation over a subset of the [L m i n , L cr i t ] inter- 
val, as in the case of the continuous cable. Courtemanche 
et al. [l(| in their analyses of a delay-integral model rep- 
resenting reentry on a 1-D loop have predicted the ex- 
istence of an infinite number of QP modes, with spatial 
wavelengths near L cr n given by 



VI. DISCUSSION AND SUMMARY 

Increasing R in the discrete loop allows sustained reen- 
try to be maintained in much shorter circuits than in 
continuous loops with equivalent lumped resistance. The 
locus of the bifurcation from period- 1 to QP propagation 
can still be predicted from the A(D) dispersion curve 
constructed by gathering data from P-l solutions and 
from mode-0 QP solutions close to the supercritical bi- 
furcation. However, increasing R modifies A(D) and the 
value of Peru ■ On one hand, the latency of the cell to cell 
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propagation is increased, due to the decrease of the in- 
tercellular current. This prolongs D because it includes 
the latency and shifts A(D) to the left. On the other 
hand, once all neighboring cells are activated, they ex- 
change current that decreases the differences in potential 
induced by the latencies. That prolongs A and, as a con- 
sequence, lifts the A(D) curve. Both effects act together 
to increase P cr it = Dcrit + A(D cr i t ). In order to analyze 
these effects, it would be more appropriate to separate 
the latency from the diastolic interval, redefining D from 
the end of the action potential to the minimum of V in re- 
polarization and considering the latency Lot from the end 
of D to the maximum derivative of the upstroke. Then 
both A and Lot could be analyzed as functions of D and 
R. However, even with this change, it will be difficult to 
build a low dimensional equivalent model of the propa- 
gation, extending the integral-delay model developed for 
the continuous loop. As seen in Fig [5] increasing R en- 
hances the dual structure of A(D) during propagation, 
which results from the effect of coupling in repolarization. 
Moreover, the same type of dual structure also appears 
for 1/9, which is almost equivalent to the latency at 
high R values. It thus becomes impossible to neglect the 
modulating effect of coupling on both A and 9 at high R 
values. Whether alternative approaches that have been 
proposed for the continuous loop would be more appro- 
priate remains to be determined. [Tl|, [Hi Gil I n an y case, 
we are still far from a general low-dimensional model that 
could be applied in situations inclu ding a dynamic change 
of the intercellular coupling, as in [281 129. l3d|. 

R also influences L m i n , the minimal length with QP 
propagation, because higher R necessitates more ex- 



citability for propagation. As a consequence, the min- 
imum D allowing sustained QP reentry increases with R 
until it reaches D cr i t (R). From this value of R, QP prop- 
agation becomes impossible. For R above this limiting 
value, period-1 reentry ends abruptly when its D reaches 
the minimal value allowing propagation. The minimal 
L for propagation increases until reaching the value of R 
where propagation becomes impossible even in a medium 
at rest. Again, the dynamics would be very interesting 
to study in a medium with dynamic modulation of the 
gap resistance. 

In all cases with QP propagation, we found the bifur- 
cation from period-1 to mode-0 propagation to be super- 
critical. For some R values, the nature of the bifurcation 
was further ascertained by prolonging simulation up to 
100 rotations. As in the continuous case, the mode-1 so- 
lutions were found to exist in a subset on the [L m i n , L cr n\ 
interval. We also devoted much effort to find n > 1 modes 
of QP propagation for different values of R, building ini- 
tial conditions either from mode-0 or mode-1 solutions 
for different L within the [L m i n , L cr i t ] interval. All these 
attempts were unsuccessful. Our initial guess was that 
the increase of R should allow more abrupt gradients of 
potential to exist between the cells, thus permitting the 
existence of higher modes of propagation. However, as 
seen in Fig [5] , the dual structure of the A(D) and 1/9 
relations becomes more pronounced at high R. This sug- 
gests that coupling still limits the gradient below what 
would be needed for higher modes of propagation. 

This work was supported by a grant of the Natural 
Sciences and Engineering Research Council of Canada. 
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